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We theoretically study an interacting few-body system of Rashba spin-orbit coupled two- 
component Bose gases confined in a harmonic trapping potential. We solve the interacting Hamilto- 
nian at large Rashba coupling strengths using Exact Diagonalization scheme, and obtain the ground 
state phase diagram for a range of interatomic interactions and particle numbers. At small particle 
numbers, we observe that the bosons condense to an array of topological states with n + 1/2 quan- 
tum angular momentum vortex configurations, where n — 0, 1, 2, 3... At large particle numbers, we 
observe two distinct regimes: at weaker interaction strengths, we obtain ground states with topolog- 
ical and symmetry properties that are consistent with mean-field theory computations; at stronger 
interaction strengths, we report the emergence of strongly correlated ground states. 
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I. INTRODUCTION 

Ultracold atomic gases offer an exceptional platform 
to explore many-body quantum phenomena due to out- 
standing experimental control over interatomic interac- 
tions, system geometry, density and purity QQ . Numerous 
research groups have, for example, successfully demon- 
strated the manifestation of few-body bound states and 
superfiuid states in Bose and Fermi gases in trapped atom 
experiments [21 IS]- Furthermore, phenomenal experimen- 
tal progress has been achieved with atomic gases loaded 
in optical lattices to emulate traditionally condensed- 
matter phenomena like superfluid-insulator transition, 
anti-ferromagnetism, and frustrated many-body systems 
[3H6]. However, due to the neutral nature of atomic 
gases, most experimental systems were limited to explor- 
ing quantum phenomena that would occur in the absence 
of electromagnetic fields. Recently, even this limitation 
was overcome, when laser fields were used to successfully 
generate effective magnetic and electric fields in neutral 
atoms [7]. The introduction of (synthetic) gauge fields 
in ultracold neutral atomic systems has thus opened the 
possibility of exploring a whole new set of phenomena 
that would manifest in the presence of abelian and non- 
abelian vector potentials |S]. 

In the presence of synthetic gauge fields in trapped ul- 
tracold bosonic systems, experimental evidence for spin- 
orbit (SO) coupling with equal Rashba and Dressel- 
haus type strengths was reported in a seminal paper 
[Sj. Recently, commendable experimental progress has 
also been achieved towards simulating SO-coupling in ul- 
tracold fermionic systems |10j . a phenomenon critical to 
the simulation of certain topologically insulating states 
in condensed-matter systems [11 . In the presence of SO- 
coupling, a generic Hamiltonian may be broadly classified 
in two classes: (a) one that breaks T (time-reversal) sym- 
metry, and which can be shown to be gauge-equivalent 
to a Hamiltonian in the combined presence of abelian 



and non-abelian vector potentials. For example, authors 
in Ref. |12| consider an SO-coupling Hamiltonian in the 
presence of a real (abelian) magnetic field and attempt to 
simulate the physics of traditional quantum Hall systems; 
(b) one that preserves T symmetry, and which can be 
shown to be gauge-equivalent to a Hamiltonian in a pure 
non-abelian vector potential. In this work, we study an 
SO-coupling Hamiltonian of the latter class, and discuss 
the emergence of ground states with unique topological 
and correlation properties. 

In this manuscript, we study an interacting few-body 
system of two-component Bose gases confined in a two- 
dimensional (2D) isotropic harmonic trapping poten- 
tial with Rashba SO-coupling. The manuscript is or- 
ganized as follows: In Sec. [H] we outline the model 
Rashba SO-coupling Hamiltonian and discuss various 
symmetries. We show that the Hamiltonian is gauge- 
equivalent to particles subject to a pure non-abelian vec- 
tor potential that preserves T symmetry. Then, we con- 
sider the non-interacting limit of this Hamiltonian, and 
discuss single-particle solutions at small and large SO- 
coupling strengths. We proceed to discuss the imple- 
mentation of Exact Diagonalization (ED) scheme to ob- 
tain the low-energy eigenstates of the interacting Hamil- 
tonian in the regime of interest to us - at large SO- 
coupling strengths. Then, we introduce various analysis 
techniques, namely:- energy spectrum, density distribu- 
tion, single-particle density matrix, pair-correlation func- 
tion, reduced wavefunction, entanglement spectrum, and 
entanglement entropy. Each technique would offer its 
unique perspective to the overall understanding of the 
ground state properties. 

In Sec. |III[ we discuss the phase diagram and analyze 
the ground state properties of the interacting Hamilto- 
nian at different particle numbers N, and at varied inter- 
atomic interaction strengths. At small particle numbers 
with N — 2, we illustrate the unique topological and sym- 
metry properties of ground states. In the relatively large 
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particle number scenario with N = 8, we observe that the 
ground states fall into two distinct regimes: (a) at weak 
interaction strengths {mean-field-like regime), we observe 
ground states with topological and symmetry properties 
that are also obtained via mean-field theory computa- 
tions; (b) at intermediate to strong interaction strengths 
{strongly correlated regime), we report the emergence of 
strongly correlated ground states. We proceed to illus- 
trate the topological, symmetry and strong correlation 
properties of these ground states. Finally in Sec. |IV[ we 
summarize and present concluding remarks. 



II. THEORETICAL FRAMEWORK 

A. System under study 

We study a two-component Bose gas confined in 
a 2D isotropic harmonic trapping potential: V{p) = 
Muj\{x 2 + y 2 )/2 = Mlo\p 2 /2. We consider the Rashba 
SO-coupling term, that couples pseudo-spin- 1/2 degree 
of freedom and linear momentum, of the form: Vso — 
—i\n{o x dy — (jyd x ), where \r is the Rashba SO-coupling 
strength and & x ,y,z are 2x2 Pauli matrices. The model 
Hamiltonian for the interacting system is then given by: 
U=J dr[Ho+H in t], 
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where r = (x,y) and ^ = [9^{r), i S^{r)) T denotes the 
spinor Bose field operators. The chemical potential \i is 
to be determined by the total number of bosons N (i.e., 
J dr^fy = N). For simplicity, we have assumed that the 
intra-component interaction strengths arc equal, so that 
fftt = 9ii = 9- The Hamiltonian is invariant under sym- 
metry operations associated with the anti-unitary time- 
reversal operator T = icr y C, and the unitary parity op- 
erator V — <t z I, where C and X perform complex conju- 
gation and spatial inversion operations respectively. The 
Hamiltonian is also invariant under the combined VT 
operator, which is unitary since operators V and T anti- 
commute, i.e., since [P,T]+ = 0. We further note that 
Rashba SO-coupling term breaks inversion symmetry. 

In experiments, the two-dimensionality can be real- 
ized by imposing a strong harmonic potential V{z) — 
Muj 2 z 2 /2 along axial direction in such a way so that 
H,ksT <C hu> z . For the realistic case of 87 Rb atoms, 
the interaction strengths can be calculated from the two 
s-wave scattering lengths a ~ lOOa^ and a-t-i, using 
g = V&^{h 2 /M){a/a z ) and g n = ^{h 2 / M){a n / a z ), 
respectively. Here, a z — y / h/{Mw z ) is the characteristic 
oscillator length in z-direction, and as is the atomic Bohr 
radius. Note that throughout this work, we consider in- 
teraction strengths such that a z 3> a,a^\,. In another 



possible regime of strong interactions where a z ~ a, a^, 
one needs to include confinement-induced resonance in 
the calculation of 2D interaction strengths g and |13j . 

In harmonic traps, it is natural to use the trap units; 
that is, to take Huj± as the unit for energy, and the 
harmonic oscillator length = \Jhj (Mwj_) as the 
unit for length. This is equivalent to setting h = 
ks = M = ujj_ = 1. For the SO-coupling, we intro- 
duce an SO-coupling length a\ — h 2 /{M\n) and con- 
sequently define a dimensionless SO-coupling strength 
\so = a ±/ a A = \J (Af/S 3 )Afl/y / a;_L. In a recent ex- 
periment [S], a spinor (spin-1) Bose gas of 87 Rb atoms 
with F = 1 ground state electronic manifold is used to 
create SO-coupling, where two internal "spin" states are 
selected from this manifold and labelled as pseudo-spin- 
up and pseudo-spin-down. This gives an effective spin- 
1/2 Bose gas. In this SO-coupled spin-1/2 BEC, \ S o 
is about 10. In a typical experiment for 2D spin-1/2 
87 Rb BECs [T3] , the interatomic interaction strengths are 
about g(N-l)m g n (N-l) = 10 2 - 10 3 (^ x ai)- These 
coupling strengths, however, can be precisely tuned by 
properly choosing the parameters of the laser fields that 
lead to the harmonic confinement and the SO-coupling. 



B. Gauge-equivalent form of Ho 



A generic single-particle Hamiltonian may be written 
\P, (1) in the form H g = (p — A) 2 /2M, where p = Kk is the par- 



ticle momentum and k is the wave- vector. The vector po- 
tential A may possibly have components in both physical 
space and spin space. Depending upon the commutation 
properties of the components of A, we may hence have 
an abelian or non-abelian type vector potential. The pri- 
mary motivation behind deriving a gauge-equivalent form 
is to map our model Hamiltonian "Ho onto Tig, and hence 
derive the nature of A. It is conceivable that depend- 
ing upon the nature of Ho, A could comprise of purely 
abelian components, or purely non-abelian components, 
or a combination of both. 

In order to map Hq onto Jig, it suffices to compare 
T~L g with the terms —h 2 V 2 /2M — i\ii{a x d y — <J y d x ) in 
%q. The latter terms may actually be rewritten as 
|p| 2 /2M + Xji{kya x — k x a y ). For a two-component Bose 
gas confined in a 2D isotropic harmonic trap, we have 
a two-component vector potential A, with A x ,A y be- 
ing 2x2 matrices. Comparing "Ho with H g , we ex- 
pect A x oc a y and A y oc — a x . Specifically, it can be 
shown that the vector potential is A = {A x ,A y ,0) — 
{hMu>±) 1 / 2 \go(cty, — a x , 0). In trap units, we then sim- 
ply have A = Xso^y, —cr x , 0). The term involving |A| 2 
is a constant, and can be gauged out without loss of gen- 
erality. Therefore, the strength of the non-abelian vector 
potential proportionally determines the strength of SO- 
coupling. It is further evident that [j4 a ,A B ] ^ 0, and 
that A is a pure non-abelian vector potential. Further- 
more, the T operator commutes with the SO-coupling 
term \R(k y a x — k x a y ). In essence, the model Rashba 
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SO-coupling Hamiltonian in Eqn. (|T|) is gauge-equivalent 
to particles subject to a pure non-abelian vector potential 
that preserves T symmetry. Proposals to realize vector 
potentials of similar forms have been addressed by mul- 
tiple groups [Sl fTSHTfj . 



C. Single-particle solutions 

We solve the model Hamiltonian H in the absence of 
interatomic interactions and obtain the single-particle so- 
lutions. Rewriting the Hq component in Eqn. (fTl, the 
single-particle wavefunction </>(r) = [<jy^ (r) , (f>± (r)P with 
energy e is given by 



H-osc 


-iX R (d y + id x ) 








i^R.{d y - id x ) 








= e 



where U osc = -h 2 V 2 /{2M) + V (p) 



(3) 

In polar coordinates 



(p,<p), we have —i(d y ±id x ) = e* llp [±d/dp~ (i/ p)d/dip}. 
The single-particle wavefunction takes the form 



<^m(r) 



<Mp) 
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(4) 



with well-defined total angular momentum j z , that is 
a sum of orbital and spin angular momenta. In gen- 
eral, we may denote the energy spectrum as e nm , where 
ri = (0, 1, 2...) is the quantum number for the transverse 
(radial) direction. 

The single-particle wavefunction <p m (r) is an eigenstate 
of the unitary V operator: 



<f>t(p) 

-h(pY v 



2tt 



(-l) m <Mr). 



The T symmetry preserved by the Hamiltonian results 
in a two-fold degeneracy (Kramer doublet) of the en- 
ergy spectrum: any eigenstate </>(r) — [<p*(r), </>^(r)] T 
is degenerate with its time-reversal partner l~(f>(r) — 
[</>|(r), — (^|(r)] T . This symmetry is preserved even in 
the presence of interatomic interactions, as the terms in 
interacting Hamiltonian %i n t are T-invariant. The su- 
perposition state, of ^ m (r) and its time-reversal partner 
state, is an eigenstate of the unitary VT operator: 

VT[<P m (r)+T<P m (r)} = (-l) m+1 [<Mr) +T0 m (r)]. 



We solve the single-particle spectrum by adopting a 
numerical basis-expansion method, details of which are 
outlined in our earlier work |18j . In Fig. [I] we show wave- 
functions of single-particle eigenstates at representative 
values of small and large SO-coupling strengths. It is evi- 
dent that a larger SO-coupling strength leads to increased 
oscillations and increased localization at radii determined 
by |m| in the radial direction. Corresponding wavefunc- 
tions 4>i{p) a l so have similar characteristics. In Fig. [2] 
we show the energy spectrum for single-particle states at 
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Figure 1: (color online). Plots (a) and (b) show wavemnctions 
</>t(p) of single-particle states in the n — manifold at small 
and large SO-coupling strengths respectively, m — (solid 
black), m = 1 (dotted red), m = 2 (dash-dotted black) and 
m = 3 (dashed red). 
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Figure 2: (color online). Plots (o) and (b) show energy 
spectrum of single-particle states at small and large SO- 
coupling strengths respectively: n = — > 7 (bottom— >top) 
and m = —16 — > +15. While energies of states within each n 
are represented by a specific symbol, it is evident that states 
with higher n have progressively higher energies. 



small and large SO-coupling strengths. From Fig.|2|a), it 
is evident that the energy spectrum is strongly dispersive 
in m at small SO-coupling strengths, with a large overlap 
between the energies of single-particle states with differ- 
ent radial quantum number n. Qualitatively, the energy 
spectrum at small SO-coupling strengths may be under- 
stood as a weak perturbation of the harmonic oscillator 
energy levels of the two pseudo-spin components. On the 
other hand, we observe from Fig. |2j6) that the energy 
spectrum is weakly dispersive or nearly flat in m at large 
SO-coupling strengths. For the range of m shown here, 
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there is no overlap between the energies of single-particle 
states belonging to different radial quantum numbers n, 
i.e, each n manifold represents single-particle states la- 
belled by their azimuthal angular momenta to with no 
overlap with adjacent n manifolds. Furthermore, the 
harmonic trapping potential may be qualitatively under- 
stood as a weak perturbation to the energy spectrum at 
large SO-coupling strengths of the corresponding trans- 
lationally invariant system. 

The localized nature of the wavefunctions in Fig. |TJ6) 
and the weakly dispersive nature of the single-particle 
energy spectrum in Fig. [2^6) are characteristics that jus- 
tify a comparison of the single-particle basis states at 
large SO-coupling strengths with 2D Landau Level (LL) 
structures in magnetic fields. In Ref. |16| . the authors 
discuss the mapping between Ho and 2D LL Hamilto- 
nian in a rigorous fashion and generalize the terminology 
of LLs as 'topological single-particle level structures la- 
beled by angular momentum quantum numbers with flat 
or nearly flat spectra' [TB] . Making use of this general- 
ization, we term the n = manifold as the lowest LL 
structure (LLL), n = 1 manifold as the next highest LL, 
and so on. As seen in Fig. [2jfo) , the radial quantization 
generates energy gaps between adjacent LLs of the order 
of trap energy hio±, i.e., of order unity in trap units. 

To summarize, we emphasize that the generalized LLs 
discussed here are created by a truly non-abelian vec- 
tor potential, i.e., in the absence of any real (abelian) 
magnetic fields. The strength of Rashba SO-coupling 
strength, and in-turn the flatness of the single-particle 
energy spectra can be experimentally controlled by using 
laser fields. At large SO-coupling strengths, as shown 
for Xso — 20, we obtain a nearly flat single-particle en- 
ergy spectra. In a non-interacting two-component Bose 
gas, quantum statistics obviates the occurrence of corre- 
lated states in a spectra that is not perfectly flat, due to 
potential condensation of all the particles in the lowest 
energy single-particle states, identified by j z — ±0.5, of 
the LLL (n = manifold). However, in the presence 
of inter-particle interactions, nearly flat energy spectra 
is sufficiently abled to act as an interesting playground 
to allow for the emergence of strongly correlated ground 
states. We now proceed to introduce the ED scheme to 
solve the interacting Rashba SO-coupled Hamiltonian at 
large SO-coupling strengths. 



states in the basis, the dimension of Hilbert space is 
D = (N + M - l)\/N\{M - 1)!. With M = 24, for ex- 
ample, D = 300 for N = 2, and D = 7888725 for N = 8. 
The dimension of Hilbert space grows dramatically with 
system size and hence, for practical purposes, we limit 
our configuration to a finite size. We observe that the 
solution becomes essentially exact when we consider a 
sufficient number of single-particle states. To solve the 
problem at hand, it is convenient to work with the SO 
single-particle basis: 



0t«™( r ) 



E 



<M r ) 
<M r ) 



dj, (5) 



where the field operator at is related to the single-particle 
state [^tnm(r),04.„ m (r)] T - Then, Eqns. ([!} and (2} sim- 
ply become 



H = E e i a l a i + E v ijki<4oja k ai, 

i ijkl 



(6) 



where (i,j,k,l) collectively denotes (n, to), and Vijki — 
(.9/2) [Vlli + K%] + g n V^ kl with 



yrt 
* ijkl 



v ijkl — 



v ijkl 



dr^(r)^.(r)^ fc (r)0 i; (r) 
dr<^>)<^- (r)0 Tfe (r)<?^ (r) 



(7) 



We perform the ED calculation in Fock space and 
the Hamiltonian T~L can be written as a matrix of di- 
mension D 2 , naively accounting for the possibility of 
inter-coupling every Fock state |19| . It is clear from 
the single particle solutions discussed in Eqn. pi), that 
the single-particle term e;ajetj contributes only to diag- 
onal entries of the Hamiltonian matrix, while the inter- 
action term Vijkia\ojakai contributes to off-diagonal en- 
tries as well. The enumeration of off-diagonal entries can 
be enormously simplified by accounting for a symmetry 
preserved by %: conservation of total angular momentum 
J z = X)jv Oxi as readily seen from Eqn. (|6j. If an entry 
Vijki is to be nonzero, we must have to.; + rrij = + to; 
in Eqn. Q. Using only the radial wavefunction, we have 
(provided to.; + rrij = urik + to;), 



D. Interacting few-body problem 
Diagonalization scheme 



Exact 



We solve the interacting Rashba SO-coupled Hamil- 
tonian H in Eqns. (JT|) and Q within the Configuration 
Interaction alias Exact Diagonalization scheme. In this 
scheme, we expand the interacting many-body Hamil- 
tonian in an appropriate single-particle basis (configu- 
ration) to obtain the solution. The solution becomes 
exact when we consider an infinite number of single- 
particle states. With N bosons and M single-particle 



v ijki 



v ijkl ~ 
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pdp(f> t i(p)<pfj(p)<t) t k(p)<l>ti(p) 
pdp fa (p) fa (p)(f>ik(p)(f>ii(p) (8) 
pdp(j) tl (p)(l) lj (p)(f> tk (p)(f> l i{p). 



This enables one to visualize the Hamiltonian in block- 
diagonal form, i.e., each block is a manifold compris- 
ing of Fock states with a fixed J z . Hence, the term 
Vijkialajai-ai can only couple states within the same 
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manifold, therefore resulting in a sparse Hamiltonian ma- 
trix. We solve this sparse matrix to identify the low en- 
ergy states of the system. 

As discussed in Sec. |II A| the Hamiltonian T~L preserves 
T symmetry. In a certain LL, the energies of states 
labelled j z and —j z are equal and hence, we need to 
consider both positive and negative angular momentum 
states in the single-particle configuration. This has two 
major implications: (a) computational intensity increases 
tremendously, and (b) a given configuration would never 
be sufficient to obtain a complete J z manifold, where 
all contributing single-particle states are included. We 
note here that the latter issue does not arise when the 
Hamiltonian breaks T symmetry, as in studies of rotat- 
ing trapped gases or gases subject to real magnetic fields 
[2Ull21| . In these studies, it was sufficient to consider only 
positive j z states and hence obtain complete J z mani- 
folds. In the limit of large SO-coupling strengths, if the 
interaction strengths are such that the energy contribu- 
tion from i/jnt is less than unity (in trap units), we may 
restrict ourselves to the lowest n = manifold. Within 
this LLL approximation, we may consider a sufficient 
number of single-particle eigenstates to obtain essentially 
exact low energy eigenstates. 



E. Analysis techniques 

ED scheme enables us to solve the Rashba SO-coupled 
Hamiltonian % and obtain the ground state phase dia- 
gram at various interaction strengths and particle num- 
bers. The ground states have interesting topological, 
symmetry and strong correlation properties. Here, we 
outline the details of various techniques that we use to 
analyze these properties. 



1. Energy spectrum 

First step in our analysis is to identify the total angu- 
lar momentum manifold J z to which the ground state be- 
longs. As discussed earlier, the Hamiltonian matrix has 
a block-diagonal form, with each block identified by its 
unique J z value. It is evident that each of these blocks 
can essentially be diagonalized independently. The en- 
ergy spectrum comprises of energy eigenvalues from each 
block, and the lowest eigenvalue and its corresponding 
J z may be readily associated with the ground state. De- 
generacies in the energy spectrum naturally reflect the 
degeneracies in the ground state. For example, a typical 
energy spectrum plot is shown in Fig. [3] 

Dimension of Fock space in the ground state J z mani- 
fold will be much smaller when compared to the Hilbert 
space dimension D. For a given parameter set, once we 
identify the ground state J z manifold, we can extract 
the coefficients of all Fock states from the correspond- 
ing eigenvector. In essence, we may then represent the 
ground state wavefunction as a sum of all contributing 



Fock states: = YlpLi a p^pi where rid is the dimen- 
sion of ground state J z manifold and a p is the coefficient 



of the Fock state <& p . As discussed in Sec. |II A| the in- 



teracting Hamiltonian H is invariant under two unitary 
symmetry operations, V and VT ■ With the knowledge 
of ground state wavefunction ~^Gi we are now equipped 
to determine if the ground state is an eigenstate of V or 
VT operator. 



2. Density distribution and single-particle density matrix 

With the knowledge of ^c, we are equipped to extract 
various properties of the ground state. We derive density 
distribution from the expectation value of single-particle 
density operator, written in second-quantized form as 



i(r') | <J(r-r') | ^(r)}^, (9) 



where |<fo(r)) is the single-particle state identified by in- 
dex j z in the LLL [21] . In our case, we also have an 
additional index to denote up- and down- spin compo- 
nents. Since J z is a good quantum number, the opera- 
tor a\a,j selects only one single-particle state within LLL 
approximation. As a consequence, it does not contain 
information about products of different amplitudes and 
loses information about interference pattern [21]. Hence, 
the density distribution solely preserves the information 
on individual densities: 



n(r) = (* G | p(r) | tf G ) 



M 

£ 

i=l 



&(r) | 2 O, 



(10) 



where Oi is the total ground state occupation of the 
single-particle state |0i(r)) [21 . Within the LLL ap- 
proximation, Oi are essentially eigenvalues of the diag- 
onal single-particle density matrix. Since single-particle 
states in Eqn. Q are eigenstates of V operator, it is ev- 
ident that the density distributions n(r) would be cylin- 
drically symmetric. For example, representative plots of 
Oi as a function of j z , and plots of density distributions 
are shown in Figs. [6] and [9] 



3. Pair- correlation function 

Pair-correlation functions help us analyze the inter- 
nal structure of the ground states. We write the 
pair-correlation operator (not normalized) in second- 
quantized form |21| . 

/°( r , r o) = ^0*( r )0j( r o)0fc(r)^(r o )a|a]a i aA ; . (11) 

ijkl 

In our case, we also have an additional index to denote 
up- and down-spin components. For instance, we may 
compute pair-correlation functions that determine the 
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conditional probability to find an up-spin or a down-spin, 
when an up-spin component is assumed to be present at 
a fixed point r , i.e., (n^(r )n^(r)) or (n-j.(r )r^(r)) re- 
spectively. We may choose ro to be away from the origin, 
but with a substantial amplitude of n(r). Due to angu- 
lar momentum conservation, the condition i + j = k + / 
must further be fulfilled. Computing the expectation 
value of /5(r,ro) with respect to ^g, we obtain the pair- 
correlation function as 

p(r,r ) =53X1 "p a P'^*(O^K r o)^( r )^( r o) 

ijkl pp' 

(% | a\a] ai a k \ %,). (12) 

When the wavefunction $c is an eigenstate of VT oper- 
ator, pair-correlation function illustrate the ground state 
symmetry properties. Furthermore, they reveal the cor- 
relations between up- and down-spin components in real- 
space. Pair-correlation functions at representative inter- 
action strengths are shown in Figs. [6] and [9] 



4- Reduced wavefunction 

We shall now discuss techniques to analyze if the 
ground states possess vortex structures with distinct 
topological properties. One identifying property is the 
presence of quantized values of skyrmion number, as dis- 
cussed in our earlier work |18| . However, this requires the 
computation of ground state wavefunction in real-space, 
a computationally prohibitive task for the bosonic few- 
particle system under study. Here, we discuss a viable 
approach to identify the topological nature of the ground 
state by computing the reduced wavefunction [22 : 



^rwf(r) 



*(rl, 



)' 



(13) 



Reduced wavefunction tp IW {(r) is computed with respect 
to one particle, here particle with index 1, while the re- 
maining N — 1 particles are placed at their most probable 
locations r? [22]. In our case, we also have an additional 
index to denote up- and down-spin components. With 
"0rwf.t( r ) an d V'rwf.^i") known, we can now extract phase 
information and compute a distinct topological quantity, 
vorticity, i.e., the number of phase slips from +ir to — tt 
along a closed contour. An integer- valued vorticity is an 
unambiguous way of establishing that the ground state 
is topological in nature with a distinct vortex structure. 
For example, typical phase plots revealing different vor- 
ticities are shown in Figs. [6] and [9] 



5. Entanglement measures 

We compute entanglement measures to analyze corre- 
lation properties of various ground states. In particular, 



we intend to probe the ground state correlation prop- 
erties that specifically stem from the presence of inter- 
particle interactions. To achieve this goal, we take cues 
from seminal papers in Ref. |23| . We choose a proper 
single-particle basis comprising of the set of eigenstates 
in Eqn. Q of the single-particle Hamiltonian Hq. In such 
a single-particle basis, entanglement in the ground state, 
or any non-degenerate energy eigenstate, occurs specifi- 
cally due to the presence of interactions |23J . 

The first step in discussing any entanglement measure 
is to partition the system and compute entanglement 
properties between different subsystems. As discussed 
in Sec. |II C| similar to 2D LL orbitals, the single-particle 
eigenstates at large SO-coupling strengths are fairly lo- 
calized in nature. This warrants us to consider partition- 
ing the system in orbital space |24| . The T symmetry 
preserved by the Hamiltonian rl naturally prompts us 
to partition the orbitals into two subsystems: positive 
j z states (subsystem A) and negative j z states (subsys- 
tem B). We write the ground state wavefunction in Fock 
space as ^g = 52 p =i a p^p> where 3> p is represented as 
| n_j c n_fj- < ._i)....Hj c _in J - c ). Here, represents the oc- 
cupation number of the single-particle eigenstate j z , and 
as discussed in Sec. |IID| a finite size cut-off is made at 
a certain value j c = j ZtC for computational feasibility. 
Now, we proceed to compute the bipartite entanglement 
properties between subsystems A and B, i.e., between 
the positive and negative j z states respectively. 

Orbital entanglement spectrum:- With the knowledge 
of we compute the entries of the density matrix p 
for the ground state as 



a v a 



p<-*-p, 



(14) 



where the generic density operator is p =| ^>g)(^g I- 

Now, we compute the reduced density matrix (RDM) 
Pa by tracing out the degrees of freedom of subsystem 
B, meaning pa=Ttb P- As shown in Ref. [23], occupa- 
tion numbers act as distinguishable degrees of freedom in 
characterizing entanglement in a finite system of identi- 
cal quantum particles. Hence in our study, RDM is com- 
puted by tracing out the occupation of all the negative 
j z states from the density matrix: 

(n 1/2 ....n jc | p jc (l/2, ..,j c ) | n 1/2 ....n jc ) = (15) 
( n -jc- n -i/2n' 1/2 ..n' je | p | n- je ..n- 1 / 2 n 1 / 2 ..n jc ) 

n- Jc ..n_i/2 

The RDM pa has a block-diagonal structure, with each 
block characterized by the total angular momentum 
that corresponds only to particles in subsystem A. The 
block-diagonal structure allows us to compute all the 
eigenvalues of the RDM using full-diagonalization tech- 
niques. Orbital entanglement spectrum (OES), termed 
so because the partition is defined in orbital space, is the 
plot of entanglement pseudo-energies as a function of 
jf. Here, £j = —\npf-. with pf being the eigenvalues of 
RDM pa [2S]. It is evident that £j with smaller magni- 
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tudes maximally contribute to the ground state proper- 
ties. 

Plots of OES reveal information about the occupation 
of various Fock states in a given ground state manifold, 
and in-turn the correlation properties of the ground state. 
If various Fock states $ p in the ground state J z manifold 
have similar magnitudes of a p , it results in similar RDM 
eigenvalues of pf, and in-turn. similar magnitudes of £j. 
Thus, if an OES plot reveals that & values are degener- 
ate or nearly degenerate, this is a clear manifestation of 
the correlated nature of the ground state. On the other 
hand, if the OES plot reveals that the values of £j are 
distinctly non-degenerate, the ground state is clearly not 
correlated. For example, representative OES plots are 
shown in Figs. [6] and [9] 

Entanglement entropy:- Plots of OES reveal the whole 
spectrum of eigenvalues of the RDM and help us un- 
derstand the correlation properties of the ground state. 
However, it is sometimes useful to extract just a single 
representative quantity from the RDM [25]. Entangle- 
ment entropy (EE) is such a measure that can be readily 
obtained from the set of eigenvalues pf of the RDM pa, 
and is defined as Sa = — tr[pA hip^] = —J2iPt^ n Pt- 
A higher entropy value means that the ground state is 
more homogeneously spread in Fock space, i.e., a larger 
number of Fock states <j> p make substantial contributions 
towards the ground state. A distinct advantage of an 
EE plot is that we are able to look at entropy values for 
a whole range of interaction strengths in a single plot, 
and thereby, understand correlation properties of various 
phases. For example, representative EE plots are shown 
in Figs. [4j [5j [7j andJH] 

In summary, density distribution, eigenvalues of single- 
particle density matrix, pair-correlation function and re- 
duced wavefunction would help us identify various sym- 
metry and topological properties of the ground states. 
Computation of RDM from proper single-particle basis 
enables us to extract various entanglement measures and 
allow us to analyze correlation properties that specifically 
stem from inter-particle interactions. 



III. RESULTS AND DISCUSSION 

As discussed in Sec. |II C| in the absence of interactions, 
all particles would simply condense into the two lowest 
energy single-particle eigenstates in the LLL identified by 
quantum numbers j z = ±0.5. This is due to the weak, 
but finite, dispersion in j z present in the single-particle 
energy spectrum shown in Fig. [2^6). The "P-eigenstate, 
identified by j z — +0.5, is represented by wavefunction 
$P = [4>\(p), (f>i(p)e iv ] T /V2n. It has a half-quantum 
vortex configuration, as the spin-up component stays in 
the s-state and the spin-down component is in the p- 
state |181 l2"71 125] . The resulting spin texture of this topo- 
logical state is of skyrmion type [18]. The degenerate 
time-reversed P-eigenstate, identified by j z — —0.5 and 
represented by T§v = [4>i(p)e~' tip , — tfi^(p)] T /\Z2n, also 



has a half-quantum vortex configuration. We may as 
well construct a zero angular momentum PT-eigenstate, 
from an equal superposition of opposite angular momen- 
tum P-eigenstates: $-p7-jz=o = ($-p ± T$-p) /V2- In 
the absence of interactions, either of the "P-eigenstates 
or the superposition "PT-eigenstate are degenerate. In 
addition, any arbitrary superposition of the degenerate 
P-eigenstates, which in principle need not be a VT- 
eigenstate, will also be a degenerate ground state. 



In the presence of inter-particle interactions, the 
ground state is not anymore determined solely by the 
energy contribution of the non-interacting part of the 
Hamiltonian Ho . Depending upon the strengths of g and 
g^, the energy contribution from the interacting part 
of the Hamiltonian %; nt also plays a crucial role. This 
competition can be better understood, especially at large 
SO-coupling strengths, by analyzing the single-particle 
wavefunctions and energy. As shown in Fig. [2^6), energy 
contributions due to Ho tries to keep the particles in 
states with lower value of angular momenta j z . However, 
for repulsive interaction strengths, energy considerations 
due to Hint tries to keep the particles as far away from 
each other as possible. This in-turn means that the par- 
ticles tend to occupy states with larger value of angular 
momenta, since they have a larger localization radii as 
shown in Fig. [TTfe). In essence, the ground state of the 
interacting many-body Hamiltonian is determined by the 
competition between the Ho and Hint terms. 



The simplest scenario where the competition between 
the Ho and H ln t terms, in-turn the effect of inter-particle 
interactions, clearly manifests is in an interacting prob- 
lem with N = 2 particles. For this reason, we discuss the 
results for N — 2 particles and analyze the ground state 
properties in greater detail, before proceeding to larger 
particle numbers. We solve the interacting few-body 
Hamiltonian H at large SO-coupling strengths using ED 
scheme within LLL approximation. The computational 
intensity, especially at large interaction strengths, limits 
the feasibility of this scheme to the order of N = 8 parti- 
cles |29| . In an earlier mean-field study on homogeneous 
two-component Bose gas [3U] , it was shown that the par- 
ticles condense into either a single plane-wave state (for 
g > g-j-j.) or a density-stripe state (for g < g^). Similarly, 
in our earlier related work on trapped two-component 
Bose gas |31| . depending on the relative magnitudes of 
g and g^, we show that states with distinct topological 
and symmetry properties emerge in the mean-field phase 
diagram. Taking cues from these results, in this study, we 
solve for the ground state wavefunction at various inter- 
action strengths, however fixing the relative magnitude 
9tl/d a * 0-5 or 1.5. In this section, we present the results 
at different particle numbers N, and analyze the topolog- 
ical, symmetry and correlation properties of the ground 
states using various techniques discussed in Sec. |II E| 



8 



0.51 

+ 



0.5 



0.50065 
0.5006 



-0.5 »0.5 

(a) g = 0.001; g n /g = 0.5 



0.5008 



CN 

^ 0.51 

+ 

^ 0.5 



0.5007 



Q 



® 



® • ® 



0.5 .0.5 



(b)j = 0.001; g n /g = 1.5 



-2 -1 -0.5 0.5 1 2 

Jz/N 

Figure 3: (color online). Energy spectrum for extremely weak 
interaction strengths with Xso = 20 and N = 2. Here, each 
marker (red) represents the lowest energy eigenvalue of a spe- 
cific block diagonal with a fixed value of J z . Since energy 
eigenvalues are very close, we identify the ground state ener- 
gies by circled (black) markers and further, show the zoomed- 
in plots in the inset. 
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Figure 4: (color online). Plots of (a) ground state J z /N man- 
ifolds and (b) entanglement entropy, as a function of interac- 
tion strength g with Ago — 20, N = 2,gf±/g = 0.5. For rep- 
resentative interaction strengths denoted by circled (black) 
markers, we illustrate the ground state properties in Fig. [6] 
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As discussed in Sec. |IIE 1| we analyze the energy spec- 
trum to identify the ground state angular momentum 
manifold J z , or equivalently, J z /N. In Fig. [3]ja), we no- 
tice that the ground state belongs to J z /N = mani- 
fold. We further determine that the ground state wave- 
function "J G is an eigenstate of VT operator. On the 
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Figure 5: (color online). Plots of (a) ground state J z /N man- 
ifolds and (b) entanglement entropy, as a function of interac- 
tion strength g with Ago = 20, N = 2,g^±/g — 1.5. 



other hand, we observe from Fig. |3j^6) that the ground 
state is degenerate in J z /N = ±0.5 manifolds. In ci- 
ther scenario, in Fig. (3^6), we determine that "Fg is an 
eigenstate of V operator. It is evident that, even in the 
presence of extremely weak interaction strengths, the in- 
teracting Hamiltonian picks either a "P-eigenstate or a 
■pT-eigenstate to be the ground state. Furthermore, it 
is clear that the ground state is sensitive to the relative 
magnitudes of and g. 

Figs. |4](a), |5](a): We solve the interacting Hamil- 
tonian H at various interaction strengths and identify 
corresponding ground state manifolds J z /N in Figs.[4|a) 
and [5ja) . It is evident from the phase diagram that de- 
pending on g and g^, the ground states belong to differ- 
ent J z /N manifolds. Furthermore, we determine if the 
ground state wavefunction ^g is an eigenstate of VT op- 
erator, and thereby identify whether the state belongs to 
V or VT symmetry phase. In a broader sense, it is evi- 
dent that a ground state in VT symmetry phase belongs 
to J z /N — manifold, while ground states in various 
J z /N 7^ manifolds belong to V symmetry phase. EE 
plots in Figs. |4|&) and [5^6) reveal correlation properties 
in various phases. For pedagogical purposes, before we 
explain the features in EE plots, we first discuss the sym- 
metry, topological and correlation properties of ground 
states. 

In Fig. [6] we illustrate density distributions, eigenval- 
ues of single-particle density matrix, orbital entangle- 
ment spectrum, pair-correlation functions and reduced 
wavefunctions at representative interaction strengths 
within various J z /N manifolds of Fig.[4|a). Using a simi- 
lar line of reasoning, we may understand the properties of 
ground states in Fig. [4^6). Let us now proceed to discuss 
various plots shown in Fig. [6] 

Figs. [6](al) — > |6](a4): In this top row, we dis- 
cuss the ground state properties of the VT eigenstate 
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Figure 6: (color online). Plots in each row illustrate the ground state properties at a representative interaction strength of 
Fig. ga). In the first column (from left), we show density distributions of spin-up component n-f(p) (solid green) and of 
spin-down component n_i(p) (dashed red). In the second column, we show eigenvalues Oi of single-particle density matrix as a 
function of angular momentum j z of the single-particle states \<f>i(r)). In the third column, we show corresponding OES plots 
of entanglement pseudo-energies £j as a function of J^/N, the average angular momentum of subsystem A. In the last column, 
we show contour plots (a4) and (64) that are normalized pair-correlation functions (n-t-(ro)n^(r)), with ro denoted by a (yellow) 
marker. Phase plots (c4) and (cZ4) are derived from reduced wavefunction tp Ci i(r), which is computed by fixing one of the two 
particles at their most probable locations and their corresponding radii are indicated by (yellow) markers. The closed dashed 
(blue) contour is a guide to the eye, that allows us to count the number of phase slips. 



in J z /N = manifold at g = 0.001 of Fig. ga). As 
shown in Fig. gal), the cylindrically symmetric density 
distributions n^(p) and n±(p) overlap. Being a VT eigen- 
state, it is evident from Fig. ga2) that the positive and 
negative angular momentum states are equally occupied. 
Furthermore, the time-reversal partner states identified 
by quantum numbers j z = ±0.5 are predominantly oc- 
cupied. As expected, from the corresponding OES plot 
in Fig.ga3), we observe that the predominant contribu- 
tion to the ground state is from the entanglement pseudo- 
energy & at J^/N = +0.25. From Figs.ga2) andga3), 
it is clear that the the maximally contributing Fock state 
is $737- =| % a =— 0.5 = l,"-j z =+o.5 = l)i which explains 
the overlapping density distributions of n\-(p) and n^(p) 
in Fig. gal). In Fig. ga4), we plot the (normalized) 



pair-correlation function (n-t-(ro)ni (r)) of this VT eigen- 
state. This plot illustrated the conditional probability to 
find a down-spin, when an up-spin component is assumed 
to be at a fixed point r , and reveals the presence of corre- 
lated regions (magnitude closer to 1) and anti-correlated 
regions (magnitude closer to 0). This plot illustrates the 
correlations present between up-spin and down-spin com- 
ponents that are not revealed by the cylindrically sym- 
metric density distributions. 

Figs. [6](bl) — > [6](fo4): In this second row, we dis- 
cuss the ground state properties of the VT eigenstate in 
J z /N = manifold at g = 0.065 of Fig. go). As dis- 
cussed with reference to Fig. gal), the density distribu- 
tions n^-(p) and n^(p) overlap in Fig.gbl). It is evident 
from Fig.go2) that the time-reversal partner states iden- 
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tified by j z — ±0.5 and j z = ±1.5 are almost equally oc- 
cupied. From the corresponding OES plot in Fig. [6^£i3), 
we observe that the ground state is equally occupied by 
& at J^/N = +0.25 and +0.75. From Figs. ^b2) and 
[6^&3), it is clear that the the maximally contributing 
Fock states are $737- =| rtj i= _o.5 = l,n^ =+ o.5 = 1) and 
(fr-pj- —\ n Ja . = _x.5 = 1, ?ij i= +i.5 = 1). To illustrate the 
internal structure of this VT eigenstate and the corre- 
lations between up-spin and down-spin components, we 
show the pair-correlation function in Fig. [6]f 64) . 

Figs. [6](cl) — > |6](c4): In this third row, we discuss 
the ground state properties of the V eigenstate in J z /N — 
+2.5 manifold at g = 0.07 of Fig. Qa). While the cor- 
responding ground state is degenerate in J z /N — ±2.5 
manifolds, we restrict our discussion to J z /N = +2.5 
manifold without loss of generality. The cylindrically 
symmetric density distributions n^(p) and n±(p) are dis- 
tinct, as shown in Fig.^cl). In this V eigenstate, there 
is an inherent asymmetry in the occupation of positive 
and negative angular momentum states. This is evident 
from the plot of single-particle density matrix eigenvalues 
Oi in Fig. |6jc2). This explains the presence of distinct 
density distributions in Fig. [BJ^cl). Furthermore, we ob- 
serve a peak in the occupation of eigenstate identified by 
j z = +2.5 in Fig. |6jc2). From the corresponding OES 
plot in Fig.|6j;c3), we observe that the ground state is pre- 
dominantly occupied by £j at J^/N = 2.5. To illustrate 
the internal structure of this V eigenstate, we show the 
phase plot derived from the reduced wavefunction ip c ^(r) 
in Fig. [6^c4). To better understand this phase plot, we 
take cues from plots in Figs. |6jc2) and[6]jc3). Though 
we observe from Fig. [6](c3) that the ground state is pre- 
dominantly occupied by at J^/N = 2.5, it may be 
conceived from Fig. |6jc2) that the ground state has con- 
tributions from various Fock states, for example: $-p =| 
%*=+2.5 = 2) or <Fp =| n^=+i. 5 = 1, n^ =+3 . 5 = 1) or 
Q-p =| nj z=+ o.5 — 1; ?ij\=+4.5 = 1). From the representa- 
tion of single-particle eigenstates in Eqn. Q, it is evident 
that the net orbital angular momentum of spin-up com- 
ponent in the ground state is +2 and that of spin-down 
component is +3. Correspondingly, the phase plot of the 
down-spin component in Fig. |6jc4) reveals a vorticity of 
3. We note here that the vorticity is the number of phase 
slips from +n to — 7T, i.e., when the shadowing changes 
from white to black. For convenience, we identify this V 
eigenstate as "P3, where 3 is the vorticity of the down-spin 
component. 

Figs. [6](dl) ^|6](d4): In this last row, we discuss the 
ground state properties of the V eigenstate in J z /N — 
+3.5 manifold at g = 0.26 of Fig. |4|a). The correspond- 
ing ground state is degenerate in J z /N — ±3.5 mani- 
folds, while we restrict our discussion to J z /N = +3.5 
manifold. As expected for a V eigenstate, the density 
distributions n^(p) and n±(p) shown in Fig. j6jc?l) are dis- 
tinct. In addition to the asymmetric occupation of posi- 
tive and negative angular momentum states in Fig.[6]jd2), 
we observe a peak occupation of eigenstate identified 
by jz = +3.5. From the corresponding OES plot in 



Fig.|6J<i3), we observe that the ground state is predomi- 
nantly occupied by at Jf/N = 3.5. However, it may be 
conceived from Fig. |6jd2) that the ground state has con- 
tributions from various Fock states, for example: $-p =| 
n^=+3.s = 2) or <Fp =| n^=+i. 5 = 1, ^,=+5.5 = 1) or 
$-p =| rij z=+ 2.5 = 1, ^=+4.5 = 1). It is clear that with 
increasing inter-particle interaction strengths, the parti- 
cles distribute themselves in higher angular momentum 
manifolds. Furthermore, it is evident that the net orbital 
angular momentum of spin-up component in the ground 
state is +3 and that of spin-down component is +4. Cor- 
respondingly, the phase plot of down-spin component in 
Fig. [6^<i4) reveals a vorticity of 4. For convenience, we 
identify this V eigenstate as VA. We further note that 
the phase plots of down-spin components derived for VI 
and V2 eigenstates in Fig. [5^a) exhibit a vorticity of 1 
and 2 respectively. 

Figs.[I](6),[5](6): As noted earlier, OES preserves the 
whole spectrum of eigenvalues of the RDM, and hence al- 
lows us to extract information about the occupation of 
Fock states with different subsystem angular momenta 
Jf. With our understanding of OES plots in Figs. [6] we 
now proceed to explain various features observed in EE 
plots of Figs. [4^6) and [5^6). (i) The presence of distinctly 
different slopes suggests the presence of distinct correla- 
tion properties in ground states within various phases. 

(ii) Within each phase, EE increases monotonously with 
increasing g. As discussed in Sec. |IIE~5| this results from 
an increasingly homogeneous distribution of Fock states 
in the ground state J z /N manifold, and in-turn an in- 
creased correlation. For example, to illustrate this fea- 
ture within the VT symmetric phase in Fig. [4^6), we may 
compare OES plots in Fig. |6ja3) and[6]j63) and observe 
an increased homogeneity in distribution of Fock states. 

(iii) The presence of nearly degenerate values results 
in a reduction in the slope of EE. While this feature is 
observed at larger interaction strengths within the VT 
symmetric phase of Fig. |4|&), the OES plot in Fig. [6^63) 
helps us understand this, (iv) Transition to a V symmet- 
ric phase is marked by a sharp reduction in the value of 
EE [32 . To better understand this feature, we compare 
OES plots in Fig. [6]f £>3) and [6|c3) and observe a sharp 
reduction in homogeneity of £j values, accompanied by 
a substantial drop in the minimum value of £j. In sum- 
mary, we emphasize that the knowledge of OES helps us 
understand various features exhibited by EE plots. 

In summary, it is evident that the interacting Hamil- 
tonian picks either a P-eigenstate or a "PT-eigenstate to 
be the ground state. The ground state is sensitive to the 
relative magnitudes of and g. J z /N plots allow us 
to identify various V and VT symmetry phases in the 
interacting system. With the analysis of density distri- 
butions, single-particle density matrix and reduced wave- 
functions, we illustrate ground state symmetry and topo- 
logical properties. We assert that the bosons condense 
into an array of P-symmetric topological ground states 
that have n + 1/2 -quantum angular momentum vortex 
configuration, with n = 0,1,2,3. With the analysis of 
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Figure 7: (color online). Plots of (a) ground state J z /N man- 
ifolds and (b) entanglement entropy, as a function of interac- 
tion strength g with \so = 20, N = 8,gfi/g = 0.5. For rep- 
resentative interaction strengths denoted by circled (black) 
markers, we illustrate the ground state properties in Fig. [9] 
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Figure 8: (color online). Plots of (a) ground state J z /N man- 
ifolds and (6) entanglement entropy, as a function of inter- 
action strength g with Aso = 20, TV = 8,gt±/g = 1.5. For 
the representative interaction strength denoted by a circled 
(black) marker, we illustrate the ground state properties in 
Fig-i 



single-particle density matrix, OES and pair-correlation 
functions, we illustrate the internal structure of different 
ground states in the VT symmetry phase. We analyze 
the correlation properties of the ground states with the 
help of OES and EE plots. 



N 



The detailed analysis presented above for the relatively 
simple, but rich, scenario of N = 2 particles shall be use- 



ful when discussing results at larger particle numbers. 
Even with a small increase in particle number from N = 2 
to N = 4 (not shown), we observe the non-occurrence of 
ground states in "P4 phase. As discussed in the introduc- 
tion of Sec. |III[ this can be understood as a manifestation 
of the competition between energy contributions from Jig 
and "Hint- A higher particle number increases the proba- 
bility distribution into single-particle states with smaller 
angular momenta, when compared to larger angular mo- 
menta eigenstates. We shall now proceed to consider the 
few-body system with N = 8 particles, discuss the oc- 
currence of various phases, and analyze the ground state 
properties at representative interaction strengths using 
various techniques outlined in Sec. |IIE| 

Figs. [7](a), |8](a): We solve the interacting Hamil- 
tonian % at various interaction strengths and identify 
corresponding ground state manifolds J z /N in Figs.[7]ja) 
and[8]ja). As discussed with reference to Figs. |4ja) and 
|5ja), it is evident that depending on g and g^, the 
ground states belong to different J z /N manifolds, and 
in-turn to VT or V symmetry phases. In this relatively 
larger particle number scenario, we observe that the 
ground states fall into two distinct regimes: (a) at weak 
interaction strengths (mean-field-like regime), we observe 
ground states with topological and symmetry properties 
that are consistent with mean-field theory computations 
|31| : (b) at intermediate to strong interaction strengths 
(strongly correlated regime), we report the emergence of 
strong correlations in ground states. The strongly corre- 
lated ground states are eigenstates of VT operator, and 
we additionally identify them with the label l SC. In 
Fig. [9] we illustrate the ground state properties at repre- 
sentative interaction strengths in these two regimes. 

Mean-field-like regime:- Figs. [9](al) — > |9](a4), 
|9](bl) ^|9](64): In the top row, we illustrate the ground 
state properties of the VT eigenstate in J z /N = man- 
ifold at g — 0.001 of Fig. ffla). It is evident that the 
properties in Figs. |9jal) — >|9ja4) are qualitatively iden- 
tical to their counterparts in Figs. [6^al) — > (6^a4). In 
the second row, we discuss the ground state properties of 
the V eigenstate in J z /N = +1.5 manifold at g = 0.013 
of Fig. ^a). The corresponding ground state is degen- 
erate in J z /N = ±1.5 manifolds, while we restrict our 
discussion to J z /N = +1.5 manifold. As expected for a 
V eigenstate, the density distributions n^(p) and n^(p) 
shown in Fig. (9^61) are distinct. It is evident from the 
single-particle density matrix eigenvalues in Fig. [9^62) 
that there is a peak in the occupation of eigenstate iden- 
tified by j z = +1.5. From the corresponding OES plot 
in Fig. we observe that the ground state is pre- 

dominantly occupied by at Jf/N = 1.5. To illustrate 
the internal structure of this V eigenstate, we show the 
phase plot derived from the reduced wavefunction tp c ^(r) 
in Fig. 9^64) . It is evident from the representation in 
Eqn. M I that the orbital angular momentum of spin-up 
component in the ground state is +1 and that of spin- 
down component is +2. Correspondingly, the phase plot 
of the down-spin component shown in Fig. [9^64) exhibits 
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Figure 9: (color online). Plots in each row illustrate the ground state properties at a representative interaction strength in 
Fig.[7|a) or[8ja). In the first column (from left), we show density distributions of spin- up component n^(p) (solid green) and of 
spin-down component n_i (p) (dashed red). In the second column, we show eigenvalues Oi of single-particle density matrix as a 
function of angular momentum j z of the single-particle states \(j>i(r)). In the third column, we show corresponding OES plots 
of entanglement pseudo-energies & as a function of J^/N, the average angular momentum of subsystem A. In the last column, 
we show contour plots (o4), (c4), and (d4) that are normalized pair-correlation functions (n^(ro)nj.(r)), with r ( ) denoted by a 
(yellow) marker. Phase plot (64) is derived from reduced wavefunction tp Ct i(r), which is computed by fixing 7 of the 8 particles 
at their most probable locations and their corresponding radii are indicated by (yellow) markers. Closed dashed (blue) contour 
is a guide to the eye, that allows us to count the number of phase slips. 



a vorticity of 2, and hence we identify this V eigenstate 
as V2. 

Strongly correlated regime:- Figs.[9](cl) — >|9](c4), 
^dl) |9](<24): In the third and fourth rows, we il- 
lustrate the ground state properties of the VT eigen- 
states in the strongly correlated regime at g = 0.021 of 
Fig. Qa) and g = 0.027 of Fig. ^a) respectively. At in- 
termediate to strong interaction strengths, as shown in 
Figs, ^a) and Fig. [8ja), all the ground states in this 
regime are eigenstates of VT operator in J z /N = 
manifold. As expected, the density distributions n^-(p) 
and ni(p) overlap in Figs. |9jcl) and [9^1). We ob- 
serve that the density distributions become increasingly 
flat with increasing magnitude of interaction strengths, 
g and g^. The interaction-induced correlations present 



in the ground states are revealed by the eigenvalues of 
single-particle density matrix and OES plots. From the 
plots in Figs.[9]jc2) and 9j^c?2) , it is evident that the parti- 
cles are nearly uniformly distributed across many single- 
particle eigenstates, with an equal distribution among 
time-reversal partner states. This distribution is quali- 
tatively in the opposite limit to the corresponding plots 
in the mean-field-like regime illustrated in Figs. [9ja2) 
and[9j&2). This feature is further substantiated in the 
OES plots of Figs. ^c3) tind^dZ), where a large num- 
ber of entanglement pseudo-energies £j are degenerate 
or nearly degenerate. As discussed in Sec. |IIE5[ the 
presence of a large degeneracy in entanglement pseudo- 
energies is a clear manifestation of the strongly correlated 
nature of the ground states. We further observe that 
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with increasing interaction strengths, the minima of the 
entanglement pseudo-energies shifts to larger J^/N 
values. To illustrate the internal structure and the cor- 
relations between up-spin and down-spin components of 
these VT eigenstates, we show the pair-correlation func- 
tions in Figs.|9j;c4) and [9^4). 



With our understanding of OES plots in Figs. [9] we 
may now explain various features observed in EE plots 
that help us understand the correlation properties of the 
ground states in the mean-field-like and strongly corre- 
lated regimes. As discussed with reference to Figs. |4]J6) 
and [5^6), we observe qualitatively similar features in 
TV = 8 particle case as well. The presence of distinctly 
different slopes in Figs. [7^6) and [8j&) suggests the pres- 
ence of distinct correlation properties in different ground 
states within various phases. Within each phase, EE in- 
creases monotonously with increasing g due to the pres- 
ence of increased correlations in the ground state. For 
example, to illustrate this feature within the VT(SC) 
phase, we may compare OES plots in Figs. [9jc3) and 
[9jc?3) and observe an increased homogeneity in Fock 
states. As a side note, we observe a small region of V- 
symmetric states before the transition to strongly corre- 
lated regime. These states do not possess distinct topo- 
logical or correlation properties. Without loss of gener- 
ality, we assert that these ground states merely occupy a 
crossover region prior to the transition to strongly cor- 
related regime. 



In summary, we emphasize that the ground states in 
the weakly interacting regime illustrated in the top two 
rows of Fig. [9] are mean-field-like states. Their den- 
sity distributions, pair-correlation functions and reduced 
wavefunctions may be readily related to the results from 
mean-field theory computations discussed in our earlier 
publication |31j. Within the ED scheme, we even repro- 
duce the reversal of phase symmetry between V and VT 
eigenstates that is observed with an increasing value of g, 
but with a fixed value of g-\\./g in our earlier mean-field 
study [31 . Such a correspondence between ED results 
and mean-field theory results is anticipated only when 
the ground state is predominantly occupied by one single- 
particle eigenstate (and/or its time-reversal partner), as 
revealed in Figs. [9ja2) and 9^62) . As illustrated in the 
bottom two rows of Fig. [9] the presence of a large degen- 
eracy in entanglement pseudo-energies and the distribu- 
tion of particles across many single-particle eigenstates, 
are clear manifestations of the strongly correlated na- 
ture of the ground states. Furthermore, we observe from 
Figs. [9] that the transition from mean-field-like regime to 
a strongly correlated regime is attained with only small 
variations in the magnitudes of inter-particle interaction 
strengths. We emphasize here that the pivotal reason 
behind this feature is the presence of nearly flat single- 
particle energy spectrum at large SO-coupling strengths. 



IV. CONCLUSIONS 

We systematically study an interacting few-body sys- 
tem of two-component Bose gases with 2D isotropic 
Rashba SO-coupling in a 2D isotropic harmonic trap. We 
show that the model Hamiltonian is gauge-equivalent to 
particles subject to a T-symmetry preserving pure non- 
abelian vector potential, whose magnitude proportion- 
ally determines the strength of Rashba SO-coupling. It 
is experimentally feasible to device a scheme in which 
tunable parameters, such as laser fields, can be used to 
control the magnitude of non-abelian vector potential, 
and hence simulate large SO-coupling strengths. In this 
limit of large SO-coupling strengths, we show that the 
single-particle energy spectrum is nearly flat. In the re- 
cent past, several research groups have made proposals 
to engineer quantum systems in which interactions would 
play a dominant role and the ground states would in-turn 
be strongly correlated. For example, recent proposals 
suggest schemes that would engineer nearly flat Chern 
bands to study strongly correlated fractional quantum 
Hall states in the lattice limit [33 . Though we study few- 
body Bose gases in traps, we emphasize that the intention 
with which we have identified the existence of nearly flat 
energy spectra at large SO-coupling strengths is not too 
dissimilar from the afore-mentioned line of thought. 

In our model system with nearly flat energy spectra, 
we observe that the presence of inter-particle interac- 
tions allows for the emergence of ground states with dis- 
tinct topological, symmetry and correlation properties. 
We solve the interacting Hamiltonian in different parti- 
cle number scenarios and analyze the ground state prop- 
erties with the help of energy spectrum, single-particle 
density matrix, pair-correlation functions, reduced wave- 
functions, and entanglement measures. At small parti- 
cle numbers, we show the phase diagram in Figs. [4] and 
[5] with ground states being eigenstates of either V or 
VT operator. In Fig. [6] we illustrate the ground state 
properties at representative interaction strengths in var- 
ious phases. We further assert that the bosons condense 
to an array of topological V eigenstates with n + 1/2 
quantum angular momentum vortex configuration, with 
n = 0, 1, 2, 3,. At large particle numbers, we illustrate the 
phase diagram in Figs. [7] and [8] We observe the presence 
of two distinct regimes: (a) at weak interaction strengths 
(mcan-field-like regime), we obtain ground states with 
topological and symmetry properties that are also ob- 
tained via mean-field theory computations. We justify 
this correspondence and illustrate the ground state prop- 
erties in detail in Fig. [9j (b) at intermediate to strong 
interaction strengths (strongly correlated regime) , we re- 
port the emergence of strongly correlated ground states. 
The properties illustrated in Fig. [9] demonstrate the cor- 
related nature of the ground states. 

It is interesting to inquire if the strongly correlated 
ground states that emerge in the nearly flat energy 
spectra would eventually allow for the manifestation of 
bosonic analogues of topological insulators predicted to 
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occur in traditional condensed matter systems. We em- 
phasize that in our system of trapped bosons, quantum 
statistics makes it impossible to fill up the lowest gen- 
eralized Landau level. This results in the absence of 
'sharp boundaries', which in-turn obviates the occurrence 
of states with topological order. However, this fundamen- 
tal roadblock may be circumvented when we consider a 
system of SO-coupled bosons or fermions in specially en- 
gineered optical lattices [34J. 
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